#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<gmp.h>
#include<time.h>

//algoritmo quartic convergence(1984)
//Then pk converges quartically to π; that is,
//each iteration approximately quadruples the number of correct digits.
//The algorithm is not self-correcting;
//each iteration must be performed with the desired number of correct digits of π.


int main(void){
	int i = 500;
	time_t t1, t2;	
	FILE *fp;
		
	fp = fopen("pi_input.txt", "r");
	if(fp == NULL){
		printf("Erro na abertura do arquivo!\n");
		return -1;
	}
    time(&t1);	

	/*Cria as variaveis da gmp*/
       mpf_t Pi;		//armazena o valor de pi
       mpf_t a0, b0, p0;
       mpf_t an, bn, pn;
       mpf_t aux, aux2, aux3;    //variaveis auxiliares
		
    /*Inicializa as variaveis da gmp*/
       mpf_init(Pi);
       mpf_init(a0);
       mpf_init(b0);
       mpf_init(p0);
       mpf_init(an);
       mpf_init(bn);
       mpf_init(pn);
       mpf_init(aux);
       mpf_init(aux2);
       mpf_init(aux3);
       
	/*define a precisao*/
	mpf_set_default_prec(33219280);

	/*atribui valor 0 as variaveis gmp*/
       mpf_set_ui(Pi,0);
       mpf_set_ui(a0,0);
       mpf_set_ui(b0,0);
       mpf_set_ui(p0,0);
       mpf_set_ui(an,0);
       mpf_set_ui(bn,0);
       mpf_set_ui(pn,0);

    /*atribui valor da constante pi na variavel gmp*/
		printf("Retorno do input file = %d\n", mpf_inp_str (Pi, fp, 10)); 

	printf("\nPi = ");mpf_out_str (stdout, 10, 10, Pi); printf ("\n\n");
    //a0 = sqrt(2);
    mpf_sqrt_ui (a0, 2);
    //b0 = 0;
    //p0 = 2 + sqrt(2);
    mpf_sqrt_ui (p0, 2);
    mpf_add_ui (p0, p0, 2);

    do{
        //an = (sqrt(a0) + (1/sqrt(a0)))/2;
        mpf_sqrt(aux,a0);
        mpf_ui_div (aux2, 1, aux);
        mpf_add (aux3, aux, aux2);
        mpf_div_ui (an, aux3, 2);

        //bn = ((1+b0)*sqrt(a0))/(a0+b0);
        mpf_add_ui (aux2,b0, 1);
        mpf_mul (aux3, aux2, aux);
        mpf_add (aux, a0, b0);
        mpf_div (bn, aux3, aux);

        //pn = ((1+an)*p0*bn)/(1+bn);
        mpf_add_ui (aux, an, 1);
        mpf_mul (aux2, aux, p0);
        mpf_mul (aux3, aux2, bn);
        mpf_add_ui (aux, bn, 1);
        mpf_div (pn, aux3, aux);

        i--;
		//a0 = an;
		mpf_set (a0, an);
		//b0 = bn;
		mpf_set (b0, bn);
		//p0 = pn;
		mpf_set (p0, pn);
		
    }while(mpf_cmp (Pi, pn) <= 0 || i > 0);

	time(&t2);
	  printf("Iteracoes = %d, Tempo calculado = %ld\n", (500-i), (int) t2-t1);
	  gmp_printf("Pi %.20Ff \n",pn);    
		

    return 0;
}
